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Abstract 

Time simulation of flutter, involving large local struc- 
tural changes, is formulated with a state-space model 
that is based on a relatively small number of generalized 
coordinates. Free-free vibration modes are first calcu- 
lated for a nominal finite-element model with relatively 
large fictitious masses located at the area of structural 
changes. A low-frequency subset of these modes is then 
transformed into a set of structural modal coordinates 
with which the entire simulation is performed. These 
generalized coordinates and the associated oscillatory 
aerodynamic force coefficient matrices are used to con- 
struct an efficient time-domain, state-space model for a 
basic aeroelastic case. The time simulation can then be 
performed by simply changing the mass, stiffness and 
damping coupling terms when structural changes occur. 
It is shown that the size of the aeroelastic model re- 
quired for time simulation with large structural changes 
at a few apriori known locations is similar to that re- 
quired for direct analysis of a single structural case. 
The method is applied to the simulation of an aeroelas- 
tic wind-tunnel model. The diverging oscillations are 
followed by the activation of a tip-ballast decoupling 
mechanism that stabilizes the system but may cause 
significant transient overshoots. 

Nomenclature 


[A] 

system matrix 

[Ai],[A 2 ] 

matrix coefficients of aerodynamic approx- 
imation 

b 

reference semichord 

[B] 

damping matrix of nominal structure 

[C\ 

output matrix 

[D],[E] 

matrix coefficients of aerodynamic approx- 
imation 
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[GB] 

generalized structural damping matrix 

[GK] 

generalized structural stiffness matrix 

[GM] 

generalized structural mass matrix 

M 

identity matrix 

k 

reduced frequency, wb/V 

[K] 

stiffness matrix of nominal structure 

m 

number of aerodynamic roots 

[M] 

mass matrix of nominal structure 

W, ] 

matrix of fictitious masses added to nomi- 
nal structure 

n 

number of degrees of freedom in finite- 
element model 

n J 

number of modes with fictitious masses 
taken into account for aeroelastic model- 
ing 

9 

dynamic pressure 

[<?(<*)] 

generalized aerodynamic force coefficient 
matrix 

[<?&(«)] 

approximated generalized aerodynamic 
force matrix 

m 

diagonal aerodynamic lag matrix 

s 

Laplace variable 

s 

nondimensionalized Laplace variable, sb/V 

t 

time 

V 

true air speed 

{*} 

discrete structural displacement vector 

{*«} 

aerodynamic state vector 

{y} 

vector of outputs 

[A B] 

change in damping matrix 

[AK] 

change in stiffness matrix 

[AM] 

change in mass matrix 

A 

eigenvalues of [A] 

[C] 

diagonal modal damping coefficient matrix 

{0 

generalized displacements 

[<t>] 

matrix of vibration modes in discrete coor- 
dinates 

W 

eigenvectors of [A] 

M 

undamped eigenvectors in modal coordi- 
nates 

U> 

vibration frequency 

K] 

diagonal matrix of natural frequencies 

W„ 

natural vibration modal frequency 


1 



Subscripts 

a actual case of stiffness variation 

accl accelerations 

b basic structural modes 

d direct model 

disp displacements 

/ fictitious masses 

F flutter 

r response locations 

T transition in stiffness 

vel velocities 

Abbreviations 

AFW Active Flexible Wing 

FM ficticious mass 

GAF generalized aerodynamic force 

MS Minimum State 


Introduction 

The common approach for formulating the equations 
of motion of aeroelastic systems starts with a normal 
modes analysis of the structural model 1 . Unsteady 
aerodynamic force coefficient matrices are then calcu- 
lated at various reduced frequency values to account for 
the aerodynamic coupling between these modes while 
undergoing oscillatory motion. Classical frequency- 
domain aeroelastic analysis methods 2 use the modal 
structural properties and the tabulated aerodynamic 
matrices to calculate the flutter conditions at which the 
aeroelastic system matrix is neutrally stable. The basic 
assumption of the modal approach is that the structural 
deflections of the aeroelastic system are linear combine 
tions of a limited set of low frequency vibration modes. 

Time-domain aeroelastic modeling techniques, which 
transform the equations of motion into a state-space, 
time-invariant form, call for the approximation of the 
aerodynamic matrices by rational functions in the 
Laplace domain. The order of the resulting state-space 
model is a function of the number of selected modes, 
the number of aerodynamic approximation roots, and 
the approximation formula. The main considerations in 
constructing the model are its size (which affects the ef- 
ficiency of the subsequent analyses), its accuracy, and 
the model construction efforts. Tiffany and Adams 3 
summarized and extended the most commonly used 
aerodynamic approximation methods. Among those, 
the minimum-state method of Karpel 4 ' 6 yields the 
smallest state-space aeroelastic models per desired ac- 
curacy. Being based on an iterative nonlinear least- 
square solution, the minimum-state method requires 
larger model-construction computer time, but the ex- 
tra time is minor relative to the time savings in typical 
subsequent analyses. 


Various computational schemes such as structural 
optimization, parametric studies, the investigation of 
damage effects, and structural changes during dynamic 
response, require repeated construction of the model for 
numerous structural variations. Repeated calculation 
of the normal modes and the associated aerodynamic 
matrices every time the structure changes is often im- 
practical in these cases. A more practical approach is 
to introduce structural changes without changing the 
modal coordinates. The validity of this approach de- 
pends on the structural information contained in the 
modal coordinates and the magnitude of the structural 
changes. The number of required modes increases with 
the magnitude of the allowable structural move limits. 

Keeping the modal coordinates unchanged is even 
more important in time simulation of aeroelastic re- 
sponse during which structural changes occur. The 
occurrences of structural changes define the time seg- 
ments between the model changes. The end conditions 
of one segment are the initial conditions of the following 
one. One can argue that these transition-point condi- 
tions can be transformed to a set of new coordinates. 
But this transformation is adequate only if the new co- 
ordinates can be expressed as a linear combination of 
the old ones. However, if this is the case, there is also 
no need to change the modal coordinates in the first 
place. An aeroelastic optimization procedure 7 where a 
realistic wing structure was optimized with aeroservoe- 
lastic stability constraints, demonstrated that moder- 
ate structural changes can be accommodated without 
changing the basic 25 modal coordinates. In order to 
accommodate large structural changes, the modes that 
serve as generalized coordinates must contain signifi- 
cant distortions in the vicinity of the changes. Tak- 
ing into account more modes supports this purpose but 
may result in an excessively large aeroelastic model. 

A method for accommodating large structural 
changes at a small number of structural locations, with- 
out significantly increasing the model size has been pre- 
sented in reference 8. The procedure starts with cal- 
culating a set of low frequency vibration modes with 
the structure loaded with large fictitious masses at the 
locations of anticipated structural changes. The ficti- 
tious masses cause the vibration modes to contain the 
local deformations required for an adequate accommo- 
dation of large structural changes. The fictitious-mass 
idea has also been applied in the context of substruc- 
ture modal coupling 9,10 and eigenvalue sensitivities of 
control augmented structures 11 . 

The purpose of this work is to outline the process of 
applying the fictitious-mass approach to generate effi- 
cient fixed- coordinate time-domain aeroelastic models 
for dynamic response with structural changes, and to 
demonstrate the use of these models for time simula- 
tion of flutter during which large local stiffness changes 
occur. 
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Analytical Development 


Direct Time- Domain Analysis 

Time-domain, first-order, state-space formulation of 
aeroelastic systems starts (as in the classic second- 
order, frequency-domain case) with zero-damping nor- 
mal modes analysis of a finite-element model of the 
flight vehicle. A set of low-frequency vibration modes, 
[<j>\, is then used to calculate the generalized un- 
steady aerodynamic force (GAF) coefficient matrices, 
[Q(ik)}, for various reduced frequency values. Time- 
domain modeling requires these GAF matrices to be 
approximated by rational interpolation functions in the 
Laplace (s) domain. The aerodynamic approximations 
in this work are performed by the Minimum-State (MS) 
method 4,5 . The approximation function, constrained to 
exactly match steady aerodynamics, is 


[Q(*)] = [g(0)] + [^i]5+[^2]s 2 + 

[D}(s[I)-m' l [E}s ( 1 ) 


where s = sb/V. The user defines the m x m diag- 
onal aerodynamic lag matrix [7?] and two additional 
approximation constraints (for each term) which define 
the [A\] and [A 2 ] matrices as functions of the other 
coefficients and the tabulated data. The [ D ] and [E] 
real coefficient matrices are calculated by an iterative 
nonlinear least-square procedure that fits the tabulated 
data matrices. The data terms can be weighted in the 
least-square process according to their relative impor- 
tance. The physical weighting algorithm 5,6 used in this 
work weights the tabulated data according to open- 
loop aeroelastic response characteristics at a selected 
dynamic pressure. It has been shown in various appli- 
cations that the resulting model is applicable to open- 
and closed-loop analyses with large dynamic pressure 
variations 5,6,12,13 . 

A full description of the state-space aeroservoelastic 
model resulting from the MS approximation is given in 
reference 14. The state-space equation of motion of the 
aeroelastic system without external excitation is 


{*} = [A]{x] 


( 2 ) 


where 


{*} = 



[A] 


0 [/] 0 

-[M)~ l [K] -?[M] _1 [£)) 

0 [E] M 


[M\ 


[GM] + ^[A?] 


[A'j = [GK] + q[Qm 
[B] = [GB} + $ [A,] 

where the [GM] is the diagonal generalized mass ma- 
trix found in the normal modes analysis, and [GK] 
is the diagonal stiffness matrix whose elements are 
GKi = GMiV* t . When the structural damping char- 
acteristics are calculated from structural element prop- 
erties (which is rarely the case), the normal modes are 
still calculated with zero damping, and the resulting 
generalized damping matrix is 

[GB] = [4>} t [BM (3) 

which is a full matrix. In most cases, however, the gen- 
eralized damping matrix is assumed to be the diagonal 
matrix 

[GB] = 2[(][GM][w n ] (4) 

where the diagonal matrix [£] is either chosen using 

engineering judgment (typical values are 0.005 to 0.02) 
or measured in a ground vibration test. The dynamic 
response at specified response locations is related to the 
response of the states of Eq. (2) by 

{1/(0} = [C]{*(<)} (5) 

where \C ] is based on the modal deflections at the re- 

sponse locations, [<f> r ]- Structural displacements are cal- 
culated with 


[C diap ] = [ 4>r 0 0 ] (6) 

structural velocities are calculated with 

[C„ el ] = [ 0 ^ 0 ] (7) 

and structural accelerations are calculated with 

[Cacc,] = -farpr] -1 [ K B qD } (8) 

The state response to initial conditions can be ex- 
pressed analytically as 

{*(<)} = (9) 

where [e At ] is a diagonal matrix where the A values are 
the eigenvalues of [A], and [rp] is the matrix of asso- 
ciated eigenvectors. Most of the eigenvalues and their 
associated eigenvectors are complex. Since the complex 
ones appear in conjugate pairs, the resulting {z(f)} 
always real, and is actually calculated by real-number 
algebra. 

Structural Changes by Direct Modal Coupling 
Structural changes that occur during the dynamic re- 
sponse, due to nonlinearities or application of mechan- 
ical devices, modify all or some of the finite-element 
[A/], [K] and [B] matrices by the increments [AM], 
[A K] and [A£] respectively. With the assumption that 
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the structural displacements after the changes can still 
be expressed as linear combinations of the original vi- 
bration modes, the analysis can continue with [ GM ], 
[GK] and [GB] of Eq. (2) replaced by 


[GM a ] = [GM] + [ <j>] T [AM][<t > ] 

(10) 

[GKa] = [GK] + [4>) t [AKM 

(11) 

[GBa] = [GB] + [4>] T [AB][4>] 

(12) 


The physically-weighted aerodynamic approximation 
performed for the original model is assumed to remain 
valid for the modified model. This has been demon- 
strated in stability analysis with moderate 7 and large 8 
structural changes. 

When the structural changes occur at time t = t ai 
the state response of Eq. (9) is continued with 

{*(<)} = [V-a][c A *(‘-‘*)][V-«]- , {x(< a )} (13) 

where A a and [\l> a ] are the eigenvalues and eigenvectors 
of the modified system matrix [A]. It should be noted 
that the modal coordinates do not change when the [A] 
matrix varies. Consequently, the response coefficient 
matrix [C] of Eq. (5) remains unchanged, except for 
that of the acceleration response, Eq. (8), which is 
now based on Eqs. (10)-(12). 

The flutter characteristics of the system before and 
after the structural changes are calculated by root-locus 
analyses of the respective [A] matrices with variable 
dynamic pressure. 

The Fictitious-Mass Model 

The application of the direct modal coupling ap- 
proach with large structural changes would require the 
inclusion of a large number of vibration modes. A way 
to analyze the system stability for large local struc- 
tural variations with a minimal increase of the model 
size and without changing the modal coordinates was 
presented in reference 8. The method is extended here 
to deal with time responses. The finite-element vibra- 
tion modes and the associated GAF’s are calculated 
as above for a nominal case, but with the locations of 
anticipated structural changes loaded with large ficti- 
tious masses. A set of ny low-frequency fictitious vi- 
bration modes, [4>j], is selected for further analysis. A 
coordinate transformation is then performed to “clean 
out” the fictitious masses and to form an actual “basic” 
case whose mass and stiffness matrices may differ from 
those of the nominal case by [AA/&] and [AKj,]. The 
transformation is based on the natural frequencies [w&] 
and eigenvectors [xb] (normalized to unit generalized 
masses) associated with the equation of free undamped 
vibrations in modal coordinates 

([GMf] + [4>f] T [AMt - {£y} + 

([G^]-b[^] T [A^][^]){^} = {0} (14) 


The mode shapes and GAF’s calculated for the 
fictitious-mass (FM) finite-element model are trans- 
formed to the basic case by 

[<t>b] = [<l>f][Xb] (15) 

and 

mm = bo>} T [Qi(ik))[xb] (i6) 

The basic-case mode shapes [0^] serve as a constant 
set of structural generalized coordinates throughout the 
response analysis (as was [<j>] in the direct modal ap- 
proach). When the columns of [<fo] and the associated 
frequencies [u/&] are compared to those calculated di- 
rectly for the basic-case finite-element model, they ap- 
pear in two groups. The low-frequency group of nj 
modes is practically identical to the directly calculated 
modes. The ny — n<j modes in the high-frequency group, 
whose number and nature depend on the number and 
magnitudes of the fictitious masses, do not represent 
actual normal modes, but are synthetic modes with rel- 
atively large local distortions in the vicinity of the fic- 
titious masses. In this way, the modes contain detailed 
local structural information with a minimal increase of 
the number of modes. 

The MS approximation of Eq. (I) is now calcu- 
lated with the transformed modal data, and the state- 
space Eq. (2) and the acceleration response matrix (8) 
are constructed for the basic case with [GM] = [/], 
[GK] = [u>&] 2 and [GB] = 2[Cfc][u>&] (in the case of 
diagonal generalized damping). When the structural 
damping characteristics are calculated from structural 
element properties, [GBj] is first calculated from the 
FM finite-element model using Eq. (3) with [<f>] = [<t>;], 
and then transformed to the basic case in the same way 
the GAF matrices are transformed in Eq. (16). 

Most of the response analysis from this point on is 
performed in the way shown above for the direct modal 
approach, Eqs. (9)-(13), with [<£] = [<f>b]. An exception 
is the generalized damping matrix, [GB a ] of Eq. (12), 
when large mass and/or stiffness changes appear in Eqs. 
(10) and (11). With no damping changes, [GB a ] would 
remain diagonal but would cause, in effect, damping 
coupling between the normal modes of the modified 
system 8 . In order to yield an effectively diagonal struc- 
tural damping matrix, [GB] is replaced by 

[GBa] = 2[GM a ][ Xa ][C a ][u>a][Xa] T [GM a ] (17) 

where [uj a ] and [x a ] are the complete set of frequencies 
and eigenvectors associated with 

[GM.]{t h } + [GK M ]{{>} = {0} (18) 

where [x 0 ] is normalized to yield 

[Xa] T (GM a ]( Xa ] = [/] (19) 

The accuracy of the FM aeroelastic model (which is 
based on ny modes) can be evaluated by separate com- 
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parisons of time responses of several structural variar 
tions to those obtained directly with their own nj nat- 
ural modes. Since the FM models contain more struc- 
tural information, the structural initial conditions have 
to be defined first in terms of the coordinates of the di- 
rect models. With {£d(0)} and {£<*(0)} being the direct- 
model displacement and velocity initial conditions of a 
comparison case, the equivalent initial conditions in the 
associated FM case are 

( 0 ( 0 )} = [Xa]{um ( 20 ) 

{6(0)} = [Xa]{6(0)} (21) 

where the columns of the nj xn^ transformation matrix 
[ \a ] are the n<* lowest-frequency eigenvectors associated 
with free FM modal equation of motion (18) for this 
case. 

Results 

The numerical example consists of a mathematical 
model of the Active Flexible Wing (AFW) wind-tunnel 
model 15 ’ 16 tested at the NASA Langley Research Cen- 
ter. An external ballast is connected to the tip of the 
wing through a mechanism designed to decouple the 
pitch motion of the ballast from that of the wing when 
flutter occurs. Figure 1 shows a close-up view of the 
ballast and its attachment to the wing. Details of the 
structural and aerodynamic models are given in ref- 
erence 8. A frequency-domain flutter analysis by the 
p — k method 2 showed that, with “stiff 1 or “coupled” 
pitch connection, the Mach 0.9 antisymmetric flutter 
dynamic pressure is about 1.9 psi and the flutter fre- 
quency is about 12 Hz. When the decoupling mecha- 
nism is activated, the pitch connection stiffness is re- 
duced by 96.5% (“soft” or “decoupled”) and the flut- 
ter dynamic pressure changes to about 3.0 psi and the 
flutter frequency to 31 Hz. The decoupling device thus 
provides a drastic change in both flutter dynamic pres- 
sure and the flutter mechanism. 

Three time-domain aeroelastic models, starting from 
separate NASTRAN normal modes analyses, were con- 
structed. Two were direct “coupled” and “decou- 
pled” models and the third was an FM model. The 
generalized aerodynamic force coefficient matrices in 
each case were calculated at 14 reduced frequency val- 
ues between 0 and 1.5. Minimum-State, physically 
weighted rational function approximations of the un- 
steady aerodynamics 5 were performed with 8 lag terms, 
which yielded 8 aerodynamic augmenting states. The 
physical weightings were performed at q = 1.5 psi. 
Each model had two versions, one where the approx- 
imation was constrained to match the real and imag- 
inary parts of the aerodynamic data at k — 1.5, and 
one with the real data matching constraint replaced by 
[A 2 ] = 0 in Eq. (1), to avoid repetitive inversions of 
[M] in Eq. (2). All the models were constructed with 
diagonal modal damping values of £ = 0.01. 


The FM NASTRAN model was with the “decoupled” 
pitch connection and with a fictitious pitch inertia of 
3 lb-in-sec 2 , twice that of the tip ballast, loading the 
wing end of the connection pitch spring. A set of 14 
low frequency fictitious- mass modes, [<£/], was used to 
generate the basic “coupled” model by solving for the 
natural frequencies and eigenvectors of Eq. (14), with 
[A Mb] - 0 in our case, then transforming the data to 
the “coupled” modal coordinates, Eqs. (15) and (16), 
and then performing the aerodynamic approximation 
and the model construction. A comparison between the 
first seven natural frequencies of the “coupled” model 
obtained directly from NASTRAN, those calculated 
from the 14 fictitious-mass modes and those calculated 
by direct coupling from 14 “decoupled” modes without 
fictitious masses is given in Table 1. It is clear that the 
FM model produces accurate frequencies. The flutter- 
critical ballast pitch frequency obtained from the “de- 
coupled” model without fictitious mass is 18% higher 
than the correct one, which indicates that the direct 
coupling approach is not adequate in our case. 

The basic FM aeroelastic model was used for var- 
ious stability and time-response analyses with differ- 
ent pitch connection stiffnesses between “coupled” and 
“decoupled” conditions. The transition from the basic 
“coupled” case to another case is performed by simply 
introducing the appropriate stiffness and damping cou- 
pling terms. A comparison between flutter conditions 
calculated by the frequency domain p — k method and 
those obtained from root-locus analyses of the state- 
space models is given in Table 2. The largest devia- 
tion of state-space from p — k results is in the flutter 
dynamic pressure of the “decoupled” FM model with 
[A 2 } = 0. The contributors to this error are the large 
difference between the flutter dynamic pressure (3.102 
psi) and the physical weighting dynamic pressure (1.5 
psi), the large deviation from the weighting “coupled” 
pitch spring, and the [A?] = 0 constraint. The error 
however is still below 5%. 

A root-loci plot of the low-frequency eigenvalues of 
the system matrix [A] where the pitch connection stiff- 
ness is varied in small increments from “coupled” to 
“decoupled” values at q = 2 psi is shown in Fig. 2. 
This was accomplished by varying [A] from “coupled” 
to “decoupled” linearly. The linear variations are valid 
because [ M ] in Eqs. (2) and (8) is not affected by 
the changes. It should be mentioned that with linear 
variations, while the extreme points modal dampings 
produced by the application of Eq. (17) are effectively 
diagonal, those at intermediate points are not. How- 
ever, without having specific knowledge on the damping 
behavior of the decoupling mechanism at intermediate 
points, the linear change is as good as other possible 
variations. The variation of the maximal real part of 
the eigenvalues of Fig. 2 is shown in Fig. 3. It can 
be observed that the system stability is very sensitive 
to the connection spring, especially when it is close to 
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